Capillary filling with randomly coated walls 
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The motion of an air-fluid interface through an irregularly coated capillary is studied by analysing 
the bucas- Washburn equation with a random capillary force. The pinning probability goes from zero 
to a maximum value, as the interface slows down. Under a critical velocity, the distribution of waiting 
times r displays a power-law tail ~ r~ 2 , which corresponds to a strongly intermittent dynamics, also 
observed in experiments. We elaborate a procedure to predict quantities of experimental interest, 
such as the average interface trajectory and the distribution of pinning lengths. 

PACS numbers: 47.55.nb,68.03.Cd,47.61.Jd 
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The ever-growing technological capability of shaping- 
up new micro-devices has revived a keen interest of the 
scientific community towards the problem of a liquid- 
vapor contact line moving on solid surfaces [ij. This is 
a widely studied phenomenon, which presents subtle ef- 
fects at different length-scales, challenging hydrodynam- 
ics, thermodynamics, and non-equilibrium statistical me- 
chanics. At the same time, this issue provides a case- 
study for a whole range of industrial applications, where 
few properties of the system, e.g. surface smoothness 
or chemical coating patterns, can be tuned in order to 
achieve the desired imbibition efficiency. Other authors 
in the past have studied the evolution of the contact line 
on a heterogeneous surface, focusing on the deformation 
of the line along transversal directions [2j, y, |4j. Our 
aim here is to provide both qualitative and quantitative 
results for the case of a narrow capillary with non ho- 
mogeneous walls [5|, |6|. To this purpose, we focus on 
the dynamics of the interface midpoint only, all other de- 
tails being projected out through the introduction of a 
position-dependent capillary force. The specific source 
of irregularity is not crucial (being it wall roughness or 
random chemical coating) as long as it can be described 
in terms of a fluctuating capillary force experienced by 
the fluid-vapor interface. An important difference with 
previous works is the inclusion of all inertial and viscous 
effects, which make the problem highly non-linear, thus 
leading to non trivial scenarios even in the simple case of 
finite memory randomness. 

The Lucas- Washburn equation [7|, |8| is a credited 
model to describe th e dy namics of a fluid penetrating 
an empty capillary [9|, [lfj, whose interface midpoint z(t) 
obeys the dynamic equation: 



dt 2 



-nz^ + m, 



(1) 



where, in 2D, 



n 



12 






is the effective drag, pi the 



fluid dynamic viscosity, pi its density and H the height 
of the channel [y, [ll(. The term f(z) is the capillary 
force determined by the wettability properties of the 



surface with respect to the fluid, f(z) = c g = 
2cos(9(z))V ca pVdiff, being 7 the fluid surface tension, 
V cap = l/pi the capillary speed and V d iff = pi/(piH) 
the diffusive speed. The angle 9{z) is usually approxi- 
mated by the static contact angle, which depends on the 
free energy balance of the solid-liquid-vapor contact line 
at rest. For the case of moving fronts, the static contact 
angle should acquire dynamical corrections [3J, ll2j . How- 
ever, these corrections arc expected to play a negligible 
role as compared to contact angle fluctuations originated 
by the irregular coating. The two terms on the left-hand- 
side of Eq. |T]) stem from the time derivative of the total 
momentum of the fluid (d/dt(pizHz)), which enters the 
capillary from an infinite reservoir. The dissipative term 
on the right-hand-side is due to friction with the walls, 
which is proportional to the filled length and to the inter- 
face velocity, which is taken to coincide with the average 
value of a Poiscuillc transversal velocity profile. 

The model of random coating used here consists of 
a sequence of patches of length A, such that the cap- 
illary force is constant on each patch: f(z) — fi for 
z G [zi, Zj+i], with Zi = iA and i = 0, 1, .... We take fi to 
be a random variable with (fifj) = (fD^ij- The proba- 
bility density function (pdf) of fi, Pf(fi), is independent 
of i, i.e. the random coating is stationary. In a real cap- 
illary, the force / can take values in a bounded interval, 
being proportional to cos(0(z)): we therefore consider 
/- < fi < /+, where /_ (+) = 2cos(6»_ (+) )V r CQp V r dl// . 
Since we focus on filling experiments, such that the ini- 
tial position of the interface coincides with the capillary 
inlet, we also require that /+ > 0. 

After a general discussion of the mathematical proper- 
ties of Eq. (Tl]), for illustration purposes, we shall present 
explicit calculations for Pf(f) in the case of uniform dis- 
tribution of the capillary force, although our analysis is 
by no means restricted to this specific distribution. The 
present model is inspired to a criterion of maximum sim- 
plicity: in particular, the coating has no long-range cor- 
relations (memory is lost above a length A). Despite this 
simplicity, our model is found to exhibit a very rich phe- 



nomenology, including power-law tails in the distribution 
of waiting times, which closely evokes the stick-slip be- 
havior observed in recent experiments 13|, [l4 1 . 

Upon introducing the dimcnsionlcss variables v = 
zz/(AV A ), s = r\t and g(z) = f{z)/V^, with V A = At/, 
Eq. (fl]) can be mapped onto a 'simple' relaxation equa- 
tion: 



dv 
ds 



— = -(v-g(z)), 



(2) 



where, in the following, we shall refer to the variables v 
and g as to "momentum" and "force" , respectively. The 
change of variable / — > g defines the boundaries g_ = 
f-/V^ and g + = f+/V A , as well as a transformed pdf 
-fg(ff) = Pf{f)\df / dg\. Because of the strongly non-linear 
^-dependence of g(z), Eq. ([2|) is very hard to solve with 
the standard analytical tools of the theory of stochastic 
processes. However, this equation permits to glean useful 
information on the local interface dynamics, i.e. when 
z £ [zj,2j_|_i], so that g(z) is constant. In particular, 
the following questions naturally arise: given the front 
at position Zi, with a given momentum v, what is the 
probability for the front to advance to z i+1 ? And, what 
are the corresponding "waiting time" r and velocity v' , 
once the next location Zj+i is reached? 

For any value of g, the answers to these questions are 
exactly determined. Given that g is a random variable, 
however, r and v' also inherit a stochastic character, the 
corresponding probability distributions being denoted as 
P t (t\v,i) and Pg v (dv\v,i), where Sv = v' — V is the mo- 
mentum change upon crossing the patch of length A. 
Since Sv is the increment of v and r is the increment of 
s, de-facto, these conditional distributions represent lo- 
cal propagators for the paths v(i) and s(i), i being the 
front position. 

The momentum increment Sv in crossing the i-th 
patch, characterized by the constant force g and initial 
momentum v, is given by Sv — v' — v = (g — v) (1 — e~ r ), 
where t is determined by solving the "exit equation" , 
obtained from the conditions z(0) = Zi, z(t) = Zj+i: 



gr + (v - g) (l - e~ T ) - (1 + 1/2) = 0. 



(3) 



Equation ([3]) connects the four variables i,v,T and g. 
Solving this equation with respect to r, defines a func- 
tion r(v,g;i), whose smallest positive real solutions rep- 
resent the waiting time to go from location i to (i + 1) 
for a fixed g and v. In particular, given the two force- 
extrema _g_ and g + , it is possible to define the max- 
imum (T max ) and minimum (r m j n ) waiting times, to 
exit the i-th patch. However, since t(v, g;i) is a tran- 
scendental function, these expressions can only be ob- 
tained numerically. On the other hand, Eq. (|3|) is eas- 
ily inverted with respect to g, determining the function 

g(r, v;i) = — — ll- i V e -1 e — > which is central to our discus- 
sion. In particular, given i and v, there exists a thresh- 
old g m i n (v; i) = min g(r, v\ i) < marking the minimum 

T>0 
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FIG. 1: Minimum and maximum values of the waiting time 
t as function of v, for i = 3, with g_ = —10 and g+ — 10. 
Inset: g(r,v; i) vs. r, for i = 3 and three different values of v 
(v=2 for case C, v = 9 for case B, and v = 13 for case A). 



value of the force such that the front is guaranteed to 
reach position Zj+i. For values of the force smaller or 
equal than g m in{v,i), the front is considered "pinned". 
The actual behavior of z(t) after such a pinning event 
is a damped oscillation around a position Zj with j < i 
(likely close to i). 

At a given choice of g-,g+ and i, one identifies three 
different possible situations, illustrated in Fig.[lJ depend- 
ing on the initial momentum v. A) a "conductive phase" , 
characterized by v > v P i n (i); B) a "weak pinning phase" , 
where v takes intermediate values, v cr it{i) < v < v P i n (i); 
C) a "strong pinning phase" , characterized by low values 
of momentum, i.e. v < v cr i t (i). Here the value v p i n (i) is 
obtained by inverting the relation g m in{v p i n ; i) = <?_, and 
represents the maximum momentum such that pinning is 
possible, while v cr it(i) = (i+ 1/2) is the critical value be- 
low which the front is pinned for any non-positive force. 





A 


B 


C 


V 


V > Vp in (i) 


Vcrit{i) <v< v pin (i) 


< V < V cr it(i) 


g m in{v;i) 


Qmin _^ Q— 


g- < g m in < 


gmin — U 


T-min \u'i *•) 


T~(v,g+;i) 


T(v,g+;i) 


r(v,g+;i) 


I'm ax\1J j 1) 


r(v,g-;i) 


r(v,g min \i) 


00 


Ppin(v,i) 





prob(g < g min ) 


prob(g < 0) 



TABLE I: Parameter ranges characterizing the three regimes 
A, B and C and corresponding ranges of the waiting time r 
and pinning probability p p i n . 



Each of these three situations corresponds to a distinct 
behavior of the function g(r, v; i) (see the inset of Fig.[T]), 
which is reflected into different ranges of existence of the 
waiting time r: in the cases A) and B), r is bounded both 
from above and below, while in the case C) it is bounded 



only from below. Figure [I] shows T m i n (v; i) and T max (v] i) 
for two values of i and a choice of g- and g+ . The three 
possible shapes of g(r, v; i) govern also the pinning prob- 
ability p p in(v; i): when momentum drops below the value 
v p in(i), the interface jumps from regime A to regime B 
and the pinning probability p P i n (v;i) goes from to a 
finite value. Further decreasing v, the pinning probabil- 
ity increases. When momentum v goes below the value 
v cr it(i), the interface enters the regime C, where the pin- 
ning probability p pm (v; i) = p™£ x = -g_/(g + -gJ), is at 
its maximum. Ranges for v, g m in and r, as well as pin- 
ning probabilities, are summarized in Tablc[I] Translated 
back to physical variables, the condition for the phase C, 
v < v cr i t (i), reads, at large i, as z < Va- 

The conditional pdf of the waiting times P t (t\v, i) is 
obtained from the pdf of the force, P g {x) = V£Pf(x), 
through the formula P T (r\v,i) = P g \g(T,i,v)]J(T,v;i), 
where 



J(r,v;i) 



\(vcrit(i) -v) + e t [(t + l)v - V cr it(i)]\ 
[{r-l) + e-rf 



(4) 

is the Jacobian | £ \ ■ Note that to obtain the bulk of 
P T one does not need the solution r(v,g;i) of the tran- 
scendental Eq. (J3j) . This quantity is however needed to 
retrieve the boundaries T m in an d r max . Note also that, 
when integrating between r m i„ and T max , P t (t\v,i) is 
not normalized to 1, but to 1 — Ppi n (v;i). If the inter- 
face is in phase C, the maximum waiting time is infi- 
nite: in this case one sees that J ~ (v cr it(i) — v)t~ 2 
for t — > oo. Diverging waiting times correspond to van- 
ishing values of the force g — > + . These two observa- 
tions sum up to yield a power-law tail for the waiting 
time pdf P T ~ P g (0)(v cr it(i) — v)t~ 2 : all moments (in- 
cluding the average) are divergent for this distribution. 
Such result is even more remarkable since it does not de- 
pend on the precise pdf of the force Pf, provided that 
P/(0 + ) > 0, i.e. arbitrarily small positive values of / arc 
allowed. In Fig. [2l the pdfs of r for i = 3 and two values 
of v = 3 < v cr it(i) and v — 4 > v cr i t (i), are shown. In an 
experiment with a randomly coated wall, as soon as the 
interface velocity z drops below Va , we expect to observe 
strongly fluctuating waiting times, with possible "appar- 
ent" pinning events, i.e. situations where the interface re- 
mains stuck for very long times before starting again with 
a finite velocity. Assuming A/H < 0.1, v\ ~ 10~ 6 m 2 /s 
and H = 10 _6 m, we obtain Va < 0.1m/ s, which app ears 
to be relevant to current experimental conditions 13l.ll4|. 



We continue our discussion by considering the 
pdf of the momentum increments Pg v (x\v t i) = 
J d,TP T (T\v,i)S[x — 5v(r,v,i)]. Numerical inspection of 
the analytic properties of the moment-generating func- 
tion w(\\v,i) = J dTP r (T\v,i)e X5v( - T ' v ^ reveals that 
Psv(8v\v, i) always has finite moments, even when v drops 
below v cr it(i). The conditional average increment, re- 
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FIG. 2: Probability density function of waiting times P T vs. 
r for i = 3 and two values of v = 3 < v or it(i) and v = 
4 > v C rit(i), with g uniformly distributed in [—10, 10]. Inset: 
the average increment Sv(v,i) as a function of v, for different 
values of i (1, 10, 20, 30 and 40). The dramatic emergence of 
a long-tail for the sub-critical case v < v cr u is clearly visible. 



stricted to the ensemble of unpinned trajectories, reads 



8v(v,i) = 



*(v,i) 



dTP;(r\v,i)6v(r,v,i), (5) 



TrniniV '^) 



where P*(r\v,i) = P t {t\v,i)/{1— p P i n {v,i)). This quan- 
tity is shown in the inset of Fig. [5] for different values 
of i. It is interesting to note that, for large values of 
i, 8v no longer depends on i and goes linearly with v, 
Sv ~ (v* — v), with v* a constant. This can be explained 
by assuming that, at large i, the interface is on average 
in the strong-pinning phase C, r ~ oo, so that Sv w g—v, 
i.e. v* — f Q 9+ gPg(g). The surviving (unpinned) trajecto- 
ries, tend to cluster on a constant momentum ensemble, 
z % Zi = v*A 2 rj. 

Equation ([1]) can be numerically integrated in order to 
gather a large statistics with many realizations of the ran- 
dom coating. This may take large computational time, 
mostly due to trajectories entering regime C, which can 
spend a long time in the same coated patch. A better 
alternative is to resort to numerical calculations based 
upon our analytical expressions, which can yield average 
quantities of interest for experimental or industrial de- 
sign. Starting at i = 1 with a given initial value of the 
fluid momentum, e.g. v = 0, one can iteratively gener- 
ate a "mean" trajectory v(i) = Y^j=i ^ V (^U ~ I)'.?)- We 
observe that this yields a good estimate of the average 
of the surviving trajectories (v(i)}, which, in the original 
variables, corresponds to (z) as a function of z. Such 
observation ensures that P$ v remains reasonably peaked, 
i.e. that the mean increment provides a fair estimate of 
the interface motion. 

We can next compute the average waiting time r(z) = 
J dxP T (x\v{i), i)x, based upon the aforementioned mean 
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FIG. 3: Left) An ensemble of 1000 trajectories as obtained 
from the numerical simulations of Eq. ^Q, and the average 
trajectory z(t)/A vs. gt, (solid line in the middle), obtained 
as described in the text (see Eq. ([5])). Right) Distribution 
of the pinning lengths; comparison between simulations (data 
points) , the numerical estimate of Eq. (6) (solid line) and the 
theoretical estimate e~ ql , as given in the text (dashed line). 



value of the momentum v(i). Again, by summing up 
all average waiting times, one recovers an estimate for 
the total elapsed time s(i) = Y^j=i^{j)- A satisfactory 
agreement between i(s) and a large sets of trajectories 
obtained by numerical integration of Eq. ([T]), is shown 
in Fig. \5^l. Note that this procedure makes sense only 
as long as r(i) is well defined, i.e. until v(i) > v cr it(i). 
As one can see, in Fig. [3]A. some of the simulated tra- 
jectories drop into this low-momentum state, and very 
long waiting times are observed indeed, in the form of 
long plateaux of z(t), before the interface starts mov- 
ing again. This corresponds to a sort of stick-slip be- 
havior for the front dynamics, very similar to the one 
reported in Fig. 2c of Ref. [l3( . Typical orders of magni- 
tude under experimental conditions are: V cap ~ lOOm/s, 
Vdiff ^ lQ- 1 m/s 1 V A < 0.1m/ s, yielding g ~ 10 3 cos(0). 
In Fig.fSJ'V we have taken Sg = g+—g- — 300, correspond- 
ing to fluctuations of the contact angle 5 cos(9) ~ 0.3. 

From the average trajectory v(i), an "average" con- 
ditional pinning probability p P i n (v,i) at each position i 
can also be computed. The total probability of observing 
front pinning at position i is given by the product of the 
probability of not being pinned at all locations j < i and 
the one of being pinned at j = i, that is: 



This quantity is numerically computed and compared 
with the pdf of the pinning lengths, as obtained from di- 
rect numerical simulations of Eq. {1} (see Fig. [3)3). The 
theoretical curve underestimates the pinning probability 
during the conductive phase, however, the rest of the pdf 
is well reproduced, including the initial peak due to the 
v = starting condition. 



A simpler prediction, which does not require any nu- 
merical computation of Eq. ([6]), can be obtained from 
the observation that v(i) — ► v* for i >> 1. When this 
saturation value meets the critical line v cr i t (i) ~ i, i.e. 
when i ~ v*, the pinning probability at each new patch 
is simply given by the maximum probability p™°f ■ From 
there on, pinning is just one of the two possible outcomes 
of a Bernoulli process, implying that the survival prob- 
ability decays exponentially, P p i n {i) ~ exp(— qi) with 
q = log(l — p'pl'lf)- This exponential tail, which marks 
the "end" of the capillary filling, for a uniform Pf{f), be- 
gins at position i ~ g+/2 = cos(9 + )V cap Vdiff/V£. These 
simple expressions may offer a handy way to estimate the 
pinning length under experimental conditions. 

Summarizing, we have highlighted the non-trivial 
properties of the Lucas-Washburn equation {1} with a 
random capillary force. Our analysis unveils the pres- 
ence of a regime characterized by a broad distribution of 
waiting times P T ~ t~ 2 (if /_ < 0). The same analy- 
sis also permits to develop qualitative estimates of the 
maximal pinning length z max /A ~ cos(9 + )V cap V dl f / /V% 
as well as the slope of its exponential distribution tail 
q ~ Iog[/+/(/+ - /_)] (if /_ < 0). Finally, we have 
developed a fast numerical procedure to retrieve detailed 
information, such as a more complete estimate of the pin- 
ning length distribution or of the average space-time tra- 
jectory inside the channel, given the statistical properties 
of the surface. These results could also be exploited in 
reverse, i.e. inferring information about the wall rough- 
ness statistics by performing many filling experiments on 
different capillaries. 
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